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ON THE NUMERICAL SIMULATION OF THREE-DIMENSIONAL 


TRANSONIC FLOW WITH APPLICATION TO THE 
C-141 WING 

Harvard Lomax, Frank R. Bailey, and William F. Ballhaus 
Ames Research Center 


SUMMARY 


Results computed by a finite-difference, relaxation algorithm are presented for the supercritical 
flow (A/qq = 0.825) about the C-141 airplane wing, which has sweep, taper, and twist. Comparisons 
with both wind-tunnel and flight data indicate that computed solutions of the classical transonic 
small disturbance equation can accurately simulate high Reynolds number flows when the shock 
sweep angle is small. It is also shown that this equation poorly approximates the complete potential 
equation when embedded shock waves are swept at angles greater than about 15°. Hence, a more 
consistent small disturbance equation is derived for use in more general cases. 


INTRODUCTION 


Finite-difference calculations of inviscid, transonic flows about two-dimensional airfoils and 
axisymmetric bodies have been reported by several authors (see, e.g., refs. 1 through 5; these works 
are summarized in ref. 6). Results of these calculations agree well with experimental data under the 
following conditions: (1) the presence of wind-tunnel walls are included in the calculations; {2) vis- 
cous effects are considered. 

The second condition depends on several others. First, the Reynolds number in the experiment 
must be high enough for the boundary layer (especially near shocks) to be thin. Second, the local 
Mach number upstream of a computed shock must be low enough to preclude shock-induced separa- 
tion in the experiment (M 1 & 1.3). Third, adverse pressure gradients near the trailing-edge must be 
mild enough for the confinement of what separation there is to a small region near the trailing edge. 
Finally, because of the preceding condition, comparisons between experimental and computed 
results should be made at the same lift coefficient. Fortunately, transonic aircraft fly at high Reynolds 
numbers, and they are designed to maintain flow attachment over as much of the wing as possible. 
Hence, inviscid computational methods are applicable to a variety of practical flow situations. 

One of the major problems facing aircraft designers today is the three-dimensional transonic 
interference problem. At transonic Mach numbers, three-dimensional effects due to the presence of 
a fuselage or wing tip can influence the flow field for large distances in the spanwise direction. Such 
flow fields, especially for low and moderate aspect ratio wings, are seldom close to two dimensional 
anywhere on the wing. It is pointed out, in references 14 and 15 (where the present method is 



described in detail), that treatment of the three-dimensional transonic problem involves more than a 
trivial extension of existing two-dimensional methods. Some of the major complications are: (1) the 
circulation is a function of the spanwise coordinate and must be determined as part of the solution 
process; (2) stability can be difficult to maintain, especially in supersonic flow regions; (3) swept 
and tapered planform shapes complicate application of wing boundary conditions; and (4) shocks, 
which are nearly normal to the free-stream direction in two-dimensions, can be oblique in the 
lateral direction. Compounding the last complication is the fact that oblique shocks may be poorly 
alined with the coordinate system making them difficult to capture sharply. 

A 

The effect of shock angle on shock strength is the subject of part I of this report. A comparison 
of the complete potential equation with the classical transonic small disturbance potential equation 
shows that the latter is a valid approximation only for flows with shocks that make small angles to 
the free-stream direction. The wing study presented in the parts II, III, and IV of this report applies 
to just this case, so the classical transonic equation was used for all the calculations presented 
herein. However, part I includes the development of a new small disturbance equation for more 
general cases when the oblique shock angle is greater than, say, 15°. 

Parts II, III, and IV cover the application of the numerical method to a practical case — that of 
the C-141 airplane wing. The C-141 wing was chosen for two reasons. First, it provides a good test 
for any three-dimensional numerical simulation since it has the complicating features of sweep, 
taper, camber, and twist; and second, a variety of wind tunnel and flight data are available for com- 
parison with computed results. Part II provides a description of the wing geometry and a discussion 
of the experimental data emphasizing the influence of viscosity. The numerical method is briefly 
discussed in part III, and the results are interpreted in part IV. 


PART I. DISCUSSION OF BASIC EQUATIONS 


^ Forms of the Potential Equations 

The exact nonlinear equation for inviscid, isentropic, two-dimensional potential flow can be 
writteii in cartesian coordinates in the form 


- *x 2 yi>xx-2*x*y% + ( a2 -*y 2 ^yy = ° 


(la) 


where a is the local speed of sound, and $ is the total velocity potential. In terms of total velocity 
components, this equation becomes 


(a 2 - U 2 )U X = 2 UVU y - (a 2 - V 2 )F, 


V x= U y 


(lb) 
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The expression of these equations in terms of perturbation potential and velocity components is 
treated in standard textbooks on fluid flow (see, e.g., ref. 7, p. 204). We will consider two of these 
perturbation forms in some detail. One form is referred to as the small perturbation equation and is 
derived by neglecting all terms containing products of perturbation velocity components. This 
leads to 


[ 1 - (7 + \Woo 2 4> x ]<}> xx ~ 2Moo 


2 d> y <p xy + [l-(.y-l)M 0 


C ]<Pyy~0 


(2a) 


or 


[1 - Moo 2 -( 7 + 1 )u]u x = TMoo 2 vu y - [1 - ( 7 - \)M<Ju]v y 


(2b) 


The second form is referred to as the classical transonic equation and is derived from equations (2) 
by considering the special behavior of fluid flow when the local velocity is close to the speed of 
sound. It can be expressed as 


[1 -Moo 2 - (7+ lWoo 2 0 x ]0 xx +( t>yy ~ 0 


(3a) 


or 


[l -Moo 2 - (7+ 1)W M 2 u]u x = -v 


v x = u y 


(3b) 


In equations (2) and (3), the perturbation velocities are defined by 


U/Uoo=l+u, V/Uoo = v 


(4) 


and the perturbation potential is defined such that 


u = <t> x , V = (l> y (5) 

In the following we will analyze these equations with regard to how they model characteristic lines, 
sub and supercritical flows, and shock waves. 
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Characteristics Consistently Ordered With Respect to Perturbation Velocities 

When equations (1) through (3) are used to model supersonic flows, they can be analyzed to 
find characteristic lines; that is, lines across which velocity gradients can be discontinuous. By divid- 
ing equations (lb) through (3b) into right- and left-hand sides, we imply that the right side can be 
evaluated along a given y coordinate, and the result can be used to integrate, or advance a solution 
in the x direction. One simple way to derive the equation for a characteristic line is to rotate the 
coordinates through some angle 9 C , as in sketch (a), and ask if there is a value of d c for which an 
advance along £ is not possible. The result of such a rotation for equation (lb) can be written: 



* minant of [A] is zero. From this it is well 

known and easy to show that a characteristic 
Sketch (a) line occurs i n solutions to equations (lb) when 

det[A] = a 2 -{U cos 9 C - V sin 9 C ) 2 = 0 (7) 

Since in isentropic flow 

a 2 /^oo 2 = (l/^oo 2 )- [(7- l)/2](2w + w 2 +v 2 ) (8) 

equation (7) reduces to the form 

[(1/Moo 2 )- cos 2 9c] + [2v cos 9 C sin 6 c -u(y- 1 + 2 cos 2 d c ] 

-[(7 - 1)(“ 2 + v 2 )/2 + ( u cos 9 C - v sin 6 C ) 2 ] = 0 (9) 

where the brackets separate, consecutively, terms having zero-, first-, and second-order powers of the 
perturbation velocity components. 

By an analysis identical to the one just outlined, one can show that the characteristic lines of 
equations (2) and (3) are 

(1/Moo 2 )- cos 2 6 C + 2v cos 6 C sin 6c - w(t -1+2 cos 2 6 C ) =0 (10) 

and 
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(1/A/ oo 2 ) - cos 2 6c - «(7 + 1 ) cos 2 d c = 0 


(ID 


respectively. Notice that equations (9) and (10) are identical for both the zero and first powers of 
the perturbation velocities. In view of their derivation, this is not a surprising result and serves, at this 
point, only as a check in consistency. Clearly, however, equation (1 1 ) is not in agreement with equa- 
tion (9), even through the first order in the perturbation velocities, Since, to the first-order, oblique 
shock waves bisect characteristic lines of the same family, this result can influence our choice of 
equations when we seek to compute flow fields with oblique shocks. 

We must proceed with caution, however, since equation (1) itself is based on the isen tropic 
assumption, which precludes, in the strictest sense, its use in connecting regions between which there 
is a shock having a change in entropy. (The error is quite small for many practical cases, being pro- 
portional to the third order of the flow deflection across the shock. Equation (1) and equations 
derived from it are often used to analyze transonic flows.) Nevertheless, a straightforward way to 
examine the approximation of a Rankine-Hugoniot shock that can be expected from the use of a 
given differential equation is found in the examination of the equation’s weak solution. This is pre- 
sented in a following section after some preparatory discussion regarding critical flow and sweep 
theory. 


The Critical Velocity and the Change in Equation Type 

If we express a second-order partial differential equation in the form 

A$ xx + B$ xy + C<S> yy +/($>*, <P y , <b) = 0 (1 2) 

it is said to be elliptic or hyperbolic, depending on whether the discriminant (B 2 - 4AC) is negative 
or positive, and it changes type if the discriminant goes through zero. This occurs for equation (1) 
when a 2 (Q 2 - a 2 ) goes through zero, where \Q\ is the magnitude of the total velocity. This gives 
the classical result that equation (1) changes type when the local flow speed passes through the 
speed of sound and is hyperbolic on the supersonic side. It is conventional to refer to this speed as 
the critical speed and the value of the pressure coefficient 


C P = 


P-Po 


(p u 2 )/2 


(13) 


at which it occurs, as the critical pressure coefficient or Cp*. Figure 1 shows a plot of Cp* vs. 
that applies to solutions of equation (1). The relationship is valid under the assumption that the flow 
is isentropic along a streamline between the reference Mach number M ^ and any given point. Also 
shown in figure 1, is the first-order small perturbation approximation to the critical pressure coef- 
ficient (Cp* & -2m*). This approximation follows at once by neglecting all terms with products of 
perturbation velocities in the equation for the pressure coefficient and can be expressed as 
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(14) 


-2 u**t 




It is interesting to find the conditions under which equations (2) and (3) change type. For 
equation (2a), the discriminant is 


M\y 2 - 1)0/ -M/0/7 + 1+0 1)] +(1 -Mj )-Mj 0/ = 0 (15) 


Two roots exist for (M 00 2 <f > x ), and the one of physical interest leads to the expression 


-2 u* = 


-2 

7+ 1 


- 1 - 




L^oo 2 


i-Kt-DAt+DKI-^oc 2 ) 


(16) 


where the * again indicates the value at which the modeling differential equation changes type. Notice 
that when the term with 0/ and terms of higher order, in 0y are neglected, equation (16) reduces 
identically to equation (14), the small perturbation form of the critical pressure coefficient found 
from the full isentropic equation. It is simple to show that the expression for -2n* derived from 
equation (3a) is also identical to equation (14). 

In summary, to the lowest order in perturbation velocities, equations (1 ), (2), and (3) all change 
type according to the same relationship among u, y, and as given by equation (14). 


Simple Sweep Theory 

The underlying basis for using wings with swept leading and trailing edges is generally derived 
by considering a sheared wing of infinite span and constant section (see, e.g., ref. 8). The principal 
assumption is that the velocity component parallel to an edge is constant or that the perturbation 
velocity parallel to an edge is zero. Under these conditions a shock, if it occurs, would also have to 
be parallel to an edge, and its strength would depend only on the component of velocity normal to 
the edge. This introduces the possibility of redefining the adjective “critical” in a three-dimensional 
flow so that it is related to a sweep angle, as in Lock and Bridgewater (ref. 9, p. 150). Bagley 
(ref. 8, p. 16) has gone so far as to use the term subcritical flow to mean shock-free flow, regardlessof 
the local Mach number. In this report the word critical designates a locally sonic condition or, 
as discussed in the preceding section, the condition at which the governing equation changes type. 
However, we respect the need for some notation that encompasses the effect of wing sweep to pro- 
duce wing sections with shock-free supercritical flow. For this purpose, in studying equation (1) 
the symbol Cp* is used, as before, to designate the pressure coefficient in an isentropic flow at which 
the local velocity is sonic (the “conventional” definition); and the symbol Cp*(d ) is used to designate 
the pressure coefficient in an isentropic flow at which the component of the total velocity normal to 
some plane, fixed by 6 = constant, is sonic. In the following discussion this plane is vertical and is 
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related to 6 in the manner shown in sketch 
(b). Note that Cp*( 0) = Cp*, and any variable 
with an asterisk superscript is defined accord- 
ing to the same logic. 

Thus, if q(6) is the component of velocity 
normal to an isobar plane, classical sweep 
theory leads at once to 


[q*\8)-Uj sin 2 0] =a 2 


(17) 


for the critical value of q(0). From this one can 
easily show 


1 - 


q*m 



-2 cos 2 6 


7+1 \M* cos 2 9 


Sketch (b) 

->) 


(18) 


Equation (18) can be used to derive Cp*(8 ) without approximation and the result 

,7/7-1 


C p *(9) = 


jM 0 


, , 1 + ~~~ Mj 1 cos 2 d 
7 + 1 V 2 


- 1 


(19) 


is illustrated in figure 2 for two free-stream Mach numbers. If the small perturbation approximation 
is applied to equation (18) it simplifies at once to 


-2 m*(0)« 1 - 


'g*mV 

Qoo 


-2 cos 2 d I 1 
T + 1 \ M oj, 2 cos 2 6 


( 20 ) 


Figure 2 compares equation (19) with equation (20), which is the sweep-theory counterpart of 
equation (14). 


Shocks 

So far the discussion has been restricted to properties of genuine solutions 1 to the differential 
equations — that is, solutions that result in continuous variations in the velocity components. We 
next inspect the properties of possible weak solutions to these equations. A weak solution can be 
composed of two genuine solutions separated by a surface across which the dependent variables can 
be discontinuous. We refer to this surface as a shock, regardless of the form of the hyperbolic equa- 
tions being inve stigated. For the Eulerian equations, the jump conditions for the dependent variables 

1 The terms genuine and weak solutions are used in the sense defined by Lax (ref. 10). 
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across the shock are referred to as the Rankine-Hugoniot conditions. Different jump conditions pre- 
vail for an arbitrary set of hyperbolic partial differential equations. In general, the conditions that 
must be satisfied across a shock can be found by writing the partial differential equations in con- 
servation law form and connecting the differences in the conserved variables with the direction 
cosines of the surface across which they are discontinuous (see ref. 1 0 or 11). Thus if, 

bE/bx + bF/by + bG/bz +H= 0 (21) 


is the differential equation, 


(A E) cos + (A F) cos a 2 + (AG) cos a 3 = 0 


( 22 ) 


is the equation that must be satisfied across the surface representing the shock. Here A signifies 
the jump of the variable across the shock and cos , cos a 2 , and cos a 3 are the direction cosines of 
the shock with respect to the*, y, and z axes, respectively. The two-dimensional cases are especially 
simple since for them the shock degenerates to a line. Choose the xy plane; then the local slope of 
the line is given by 


cos a 2 /cos oq = -tan d 

where 6 is an angle having the same relation to the y axis as 0 C in sketch (a). 

The jump conditions for the isentropic shock embedded in equation (1) are derived and 
discussed in reference 1 2. The conditions lead to an implicit relationship for the velocity components 
on the two sides of the discontinuity, which can be written as 


G 1 u l ~ G 2 u 2 = (G x Vi - G 2 v 2 ) tan 6 (23) 

where 



and the subscripts 1 and 2 refer to conditions on the upstream and downstream side of the shock at 
any given point. 

The conservative form of equation (2) can be expressed as 
^ r I ^ 

— (1 -^oo 2 )« - —~r M * o 2 * 2 -~^ M oo 2 (y + 1)« 2 [v-Mo^iy- 1)mv] = 0 

( 24 ) 

bv bu _ 
bx by 
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and for its jump condition, straightforward algebra leads to the relationship 


1-2^oo 2 (7+1)(“ 


1 +M 2)j 


tan 2 Q - 2 M ao 2 (v 2 + u 2 tan 0) + 1 -M^ 2 - - M^ 2 (y + l)(«j + u 2 )=0 

(25) 


The jump condition for equation (3) is important and is derived in some detail. The conserva- 
tive form is 


3 

3x 


1 3 v 

(\-M 00 2 )u--M 00 2 ('y+\)u 2 + — -= 0 

2 J ay 

3v du _ q 
dx 3 y 


and the jump conditions first appear in the form 


(1 - MJ )A u - - (7 + 1 )Aw 2 = Av tan d 


Av = -A u tan 9 


Eliminating Av, and noting that A u 2 -u x 2 -u 2 2 = (u t + u 2 )Au = u Au gives 


(26) 


tan 2 9 = M^ 2 - 1 + ^ M^ 2 (7 + 1 )u 


(27) 


or 


(-2u, ) + (-2 m 2 ) 

2 


-2 

7 + 1 



(28) 


This is the condition that must be met across a shock according to the classical transonic theory 
expressed by equation (3). 


Shocks and Supercritical Flow 


In the terminology we are using, a simulated flow becomes supercritical when the differential 
equation changes type from elliptic to hyperbolic. To the lowest order in the perturbation velocities, 
we have shown that this occurs under the same condition for equations (1) through (3), namely, 
where 


-2 u = 


-2 u* 


-2 

7+1 



(29) 
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In terms of the definitions and analysis given above, two observations can be made about a simulated 
flow that has become supercritical: 

1. A shock can exist; but, without prescribed boundary conditions, there is nothing that tells 
us it must exist; that is, the simulated flow along a streamline can go supercritical and return to sub- 
critical with or without an intervening shock, regardless of the form of the nonlinear equations being 
used to model the flow. 

2. The critical condition — for example, equation (29) — by itself tells us nothing about the 
position or intensity of a shock if it does exist. Information about these properties appears in 
the weak solution of the differential equation and is resolved by a mechanism quite different from 
that which governs change of type. 

Now consider the jump condition for the shock embedded in the equations (24) and (26). 
Recall that we are considering a three-dimensional, transonic flow having the property that a shock, 
if it exists, will be nearly vertical, so we need to consider only the effect of 9, its lateral angle of 
sweep. It is quite possible to study shocks with significant inclination to the free stream in both 
directions; but, again, the added complexity is unnecessary for our purposes. 

When all simulated shocks are nearly perpendicular to the free-stream direction so that setting 
9 = v 2 - 0 is a good approximation, equations (24) and (26) represent the same result. This result 
applies to the standard two-dimensional studies, of course, and in the form of equation (28) it gives 
a condition that must be met by the average of the perturbation velocities across the discontinuity. 
Notice that as u x approaches u 2 , they both approach u*. That is, the intensity of a shock perpen- 
dicular to the free stream is centered about the critical perturbation velocity, which, in turn, is based 
on the condition that the equation changes type. For such cases, equations (2) and (3) give essen- 
tially the same results, and equation (3) is to be preferred because of its simplicity. 

If the boundary conditions are such that 
a vertical shock is oblique to the free stream, 
the shock position and intensity could be sig- 
nificantly different depending on whether 
equation (2) or (3) was used for the simulation. 

For example, return to equation (25) and con- 
sider the special case when v 2 - -2 u tan 9 
(sketch (c)), which exists when the component 
of the perturbation velocity parallel to the 
shock is zero. That is, the special case when 
the condition for simple sweep theory applies, 
and equation (25) reduces to 

(-2 u j ) + (-2 u 2 ) ^ -2 cos 2 6 I 1 

^ y + 1 \ M oq 2 cos 2 9 

On the other hand, the jump condition for equation (3), given by equation (28), does not depend on 
any special relation between u 2 and v 2 and is valid for all values of 9 as it stands. Therefore, the 




X 


Sketch (c) 
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difference between equations (2) and (3) in simulating oblique shocks is apparent from a comparison 
of equations (30) and (28). This difference is shown in figure 3 for M m = 0.85. Shown also is the 
exact value of Cp *'(6) for simple sweep theory given by equation (19). Once again, as the strength of 
the shock vanishes, u x approaches u 2 , and they both approach a value of w*(0) appropriate for the 
local conditions. Figure 3 shows that equation (3) is a very poor model for flows with shocks that 
are more than 30° oblique to the free stream. Equation (2), on the other hand, is acceptable for a 
wide range of 0 , limited only by the accuracy of the small perturbation approximation itself. Hence, 
for three-dimensional, transonic flows in which the presence of oblique shocks at moderate to large 
sweep angles is anticipated, the proper form of the transonic small disturbance equation is 


[1 -Mj -(7+1 )MJ (f> x ) 4 > xx - 7MJ <t> y <p xy + [ 1 - (7 - 1 Woo 2 <fi x ] <P yy + <t> zz = 0 


(31) 

From the above discussion we see that shocks may or may not appear in flow simulations when 
the differential equations become supercritical. If they do form, they appear as surfaces of disconti- 
nuities, the local intensity of which depends on the local surface inclination with respect to the free-: 
stream direction as well as other local and free-stream conditions. We now inspect the effect that 
these observations have on the interpretation of numerical simulations of transonic flow. 


PART II. DESCRIPTION OF THE AIRPLANE WING AND DISCUSSION 
OF EXPERIMENTAL RESULTS 


Several numerical simulations of three-dimensional wings have been reported (refs. 13-15). 
These studies have included lift and sweep-back but not taper or twist. Further, they were confined 
to low aspect ratio planforms or biconvex profiles. To test the full capability of a complete wing 
analysis, we chose to simulate the transonic flow about a wing with a fairly high aspect-ratio planform 
having twist, sweep, taper, and a blunt, cambered profile. The wing on the C-141 airplane has all of 
these features; in addition, some flight and wind-tunnel data are available for comparison. Previous 
studies of a panel model of the C-141 profile (ref. 14) indicated that the calculations are feasible if 
they can be carried out on an advanced computer such as an IBM 370/195 or a CDC 7600. The 
numerical procedure and results are presented in part III. 


Wing Geometry 

The C-141 wing planform is shown in figure 4. It is mounted on top of the fuselage and has a 
kink in the trailing edge about 43 percent of the span outboard from the root section. Two nacelles 
are mounted under the wing as shown. The fuselage, kink, and nacelles were not modeled in the cal- 
culations. The planform that was modeled is shown by the solid line in figure 5. The actual numbers 
are shown in figure 5 and result in the following basic parameters: 
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Sweep (L.E.) = 25.6° 

Taper ratio = 0.373 
Aspect ratio = 8.0 
Thickness ratio = 0.1 145 
Twist = 5° (linear) 

At the time the flow field was calculated, we were not aware that the C-141 profile varied in the span 
direction. For the computations we assumed that the profile was uniform and used the section 
coordinates at the span station y/b = 0.389 (given in table 1), which is the measurement station 
closest to the kink as shown in figure 5. 


Brief Review of Viscous Effects on Transonic Flow Experiments 

Before proceeding with a discussion of the experimental results, the following remarks should 
be made. It is very well known that low Reynolds number wind-tunnel tests of the pressure distribu- 
tion on a complete C-141 model do not agree with flight conditions (see, e.g., ref. 16), a discrepancy 
that is attributed mostly to differences in viscous scaling rather than wind-tunnel wall or flow 
effects. The purpose of reference 16 was to present a technique for overcoming this difficulty in 
wind-tunnel facilities by testing wing panel models with stream-lined end plates to obtain experi- 
mental data at much higher Reynolds numbers than are possible for complete models. This is a per- 
fectly acceptable procedure, and its success is indicated by the results presented in reference 16. 

One purpose of the present report is to suggest an alternate method for anticipating the loads 
and pressure distributions on wings at flight conditions. In this method, numerical simulations are 
performed that, in the absence of boundary-layer interaction, can be said to represent unseparated 
flows at an infinite Reynolds number. 

At the present time the difficulty in predicting the transonic flow that exists in actual flight 
from either wind-tunnel or computer simulations is generally attributed to a viscous-inviscid coupling 
of two types of turbulent boundary-layer separations (ref. 17). One of these, called pressure gradient 
or trailing-edge separation, is associated with trailing-edge conditions. It can occur at any speed, and 
depends strongly on Reynolds number. It is connected with the loss of lift that experimental results 
show relative to theoretical results for the same wing at the same angle of incidence if the theoretical 
results are based on the enforcement of the Kutta condition in a potential flow calculation. The 
other type of separation, called shock-induced separation, is generally considered to be insensitive 
to Reynolds number and depends primarily on the pressure jump across the surface of the shock. 

In reference 17, it is argued from experimental evidence that shock-induced separation cannot 
occur until a free interaction compression (just outside the boundary layer) turns the flow through 
6.6° . Given sonic conditions on the downstream side of the shock the lowest Mach number on the 
upstream side meeting this condition is about 1.32. It is further hypothesized that the shock will 
continue to induce local separation for upstream Mach numbers greater than about 1.32, although 
the angle of separation will continue to be around 6.6° for the higher Mach numbers. In view of 
these arguments and examinations of other data, a set of tentative hypotheses (also suggested in 
ref. 17) has been adopted: 
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1 . As the free-stream Mach number in a transonic two-dimensional flow is increased, a shock 
will (in general) form and commence to move downstream. The Mach number just ahead of the 
shock can increase to about 1.2 without danger of shock-induced separation. 

2. As the free-stream Mach number is further increased, the shock will move farther downstream 

and its upstream Mach number may go as high as about 1.4 without shock-induced separation, but 
this separation might occur anywhere between 1.2 £ 1 .4. The data in figure 6(b), for example, 

show a break at M » 1.3, which is presumed to be due to shock-induced separation occurring at Mach 
numbers higher than this. 

3. When the shock does separate the boundary layer, the interaction causes the shock to weaken 
and move upstream as much as 20 percent of the chord from its unseparated location. 

Figure 7 relates the local pressure coefficient to the free-stream Mach number given by the equation 


i 


C P ~ 


yM a 


[i + [(t- n/2]^ 

+ [(7-0/2 ]M 2 




for the range of conditions covered in the three hypotheses just outlined. 

In reference 17, it is shown that the two types of separation (trailing edge and shock induced) 
can amplify one another when they occur together. There appears to be good reason to believe that 
all variations of these viscous-inviscid interactions appear in the data presented below. 


The Experimental Data 

Flight and wind-tunnel data for the C-141 airplane and its model are shown in figures 8(a) and 
(b). The data are taken directly from reference 16, and represent pressure distributions on the upper 
surfaces of the wings. To inspect these data in light of the various possible separation phenomena dis- 
cussed above, we first consider the flight data for M ^ = 0.825, Re = 36X 10 6 shown in figure 8(a). 
It- is expected that some trailing-edge separation exists since this is generally the case for all wings and 
wing sections in flight or wind tunnels (typically a loss of 10 percent in lift under the conditions 
reported here). However, there is no indication of shock-induced separation. The pressure coefficient 
just ahead of the shock is about -1.07 which, according to figure 7, is just within the allowable 
limits of unseparated flow according to the second hypothesis given in the preceding section. (The 
complicating factor of sweepback also enters the picture. However, on the basis of the computed 
results the shock is probably swept no more than 13°.) 

When the flight Mach number is increased to 0.85, the local Mach number ahead of the shock 
is pushed above the “allowable” limit. From the appearance of the data, shock-induced separation 
occurs causing the shock intensity to decrease and the shock location to be moved upstream from 
its position in an inviscid flow (about 20 percent of the chord according to calculations presented in 
the next part). On the other hand, when the Mach number is fixed at 0.825 and the Reynolds num- 
ber is reduced from 36X1 0 6 to 8.5X10 6 , the data from the model show that the shock also moves 


13 



upstream. Presumably, this is attributable to trailing-edge separation, which is strongly influenced 
by Reynolds number. 

Finally, when the Mach number for the model tests is increased to 0.85, the data shown in 
figure 8(b) can be interpreted to show strongly coupled effects of both types of separation. 


PART III. DISCUSSION OF THE NUMERICAL PROCEDURE 
AND PARAMETERS 


The procedure used in the calculations presented in this report is as follows: 

1. Test the local flow condition using central differencing. 

2. If, on the basis of this test, the flow is supersonic, use upstream differencing for the stream- 
wise derivative. 

3. If, on the basis of this test, the flow is subsonic, use central differencing. 

This procedure consistently underestimates the strength of the shock and appears to cause its loca- 
tion to fall upstream of its position as determined by the exact solution of the governing partial dif- 
ferential equation. This can easily be shown to be true in the case of a one-dimensional flow equation, 
and is implied to be true by experience with two- or three-dimensional cases. 

On the other hand, the viscous effects brought about by shock-boundary-layer interaction 
also cause the experimental strength and location of the shock to differ from those values predicted 
on the basis of a theory that neglects viscosity. If the predominant viscous effect is limited to trailing- 
edge separation, then it appears by coincidence that computed flows (using methods based essen- 
tially on the procedure described above) and experimental flows are in good agreement. A substan- 
tiation of this remark can be found by comparing the computed results for two-dimensional flows 
shown in figure 6(a) with the high Reynolds number experimental results shown in figure 6(b). The 
usefulness of this fortuitous circumstance cannot be overlooked now that it has been established by 
such a wide variety of comparisons between calculated and experimental results. It is hardly a 
satisfactory state of affairs, however, and further work is required to develop methods that account 
for viscosity in a rational manner. 

Further details of the numerical method used to simulate the flow are described in references 14 
and 15. On the basis of experience gained in these calculations, the authors are developing methods 
having greater efficiency and reliability than that given in references 14 and 15. 

Part of the numerical technique was to use a coordinate transformation that transformed the 
wing into a rectangular area for the final mesh that was used in the differencing. Both a coarse 
(~80,000 points) and a moderately fine (~200,000 points) mesh were used for the calculations; 
the results are described in the next part. In both cases the chord distribution of points was 
clustered near the leading and trailing edge, and nearly equally spaced over the remaining portion 
of the wing. As the mesh receded from the wing, the intervals between points were successively 
increased in all directions. The actual array of (x, y, z ) points used was (68, 23,49) for the coarse 
mesh and (82, 49, 49) for the finer one. Other details of the calculations are: 
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Coarse mesh 

0.0025 

0.05 

-0.06 

1.3 

Fine mesh 

.0025 
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PART IV. INTERPRETATION OF NUMERICAL SIMULATIONS AND 
COMPARISON WITH EXPERIMENTAL RESULTS 


All calculations presented were based on the equation 


[1 -Mqo 2 ' (7+ \)MJ<1> X U XX +<t> yy +<p zz = 0 


(32) 


The Kutta condition was applied at the trailing edge by requiring that the pressure be continuous 
across the wake. The vortex sheet was fixed in the plane of the wing, and its spanwise strength was 
proportional to the span loading. The boundary conditions were applied by specifying the upper 
and lower vertical velocities above and below a horizontal plane according to the local thickness, 
camber, and angle of attack of the wing surface (thin-airfoil boundary conditions; see ref. 14 for 
details). In studies of two-dimensional profiles, this approach has been found to give results that com- 
pare reasonably well with data for airfoil nose shapes similar to that found on the C— 141 section 
(those having “nonpeaky” pressure distributions, see ref. 18). One would expect the approach to 
give comparable results in three-dimensional studies having twist, taper, and sweep, except for cases 
in which oblique shocks form on surfaces at an angle greater than about 1 5° with the free stream. 
The latter restriction is brought about by the use of equation (32), rather than the boundary 
condition approximation, and follows from the discussion given in part I of this report. 

It is of fundamental interest to study the numerical resolution of shocks that are admitted by 
equation (32), whether or not they represent a valid approximation of normal or oblique Rankine- 
Hugonoit shocks. As noted in part I, the interpretation of a shock is related entirely to an analysis 
of the set of differential equations and boundary conditions from which it was derived and is valid 
and applicable to that set irrespective of how that set models a physical process. Therefore, in spite 
of the fact that some of the following numerical results for the M ^ = 0.853 case may have no quan- 
titative physical counterpart, their qualitative analysis is quite instructive. 


Examples of Normal and Oblique Shocks, M ^ - 0.853 

The calculations made for M = 0.853 were carried out on a relatively coarse mesh (80,000 
points) and are presented in figures 9, 10, and 11. The computed shock location at the span loca- 
tion y/b = 0.414 occurred at x/c « 0.78. This fact, taken in conjunction with the results shown in 
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figure 8, is the basis for our conclusion that the flight test results for M = 0.853 are strongly 
affected by shock-induced separation. However, they are still worth considering from an instruc- 
tional point-of-view. 

Figure 9 presents the computed results for the pressure on the upper and lower surface of the 
wing at the root section. Notice that the curve for the pressure on the upper surface shows four 
places where there are sudden changes in curvature. To properly interpret the solution, one must 
decide whether these are valid representations of an analytic solution to equation (32) or invalid 
representations caused by an inadequate numerical approximation. Then, if they are considered to 
be valid, one must decide whether they represent an ordinary and continuous response to local 
wing shape, or a weak solution with a shock. 

Seeking out the possibility of a shock at the root section, we note that (1) a shock can occur 
at any position where Cp (= -2 u) is greater than 0.312 (where eq. (32) changes type); (2) because of 
symmetry, at this station it cannot be oblique to the y axis; and (3) because of (1) and (2), the 
average Cp across any discontinuity must also be 0.312. With this in mind, we find that figure 9 
shows two locations where the results could be associated with a shock, one at x/c = 0.1 70, and the 
other at 0.945. Before considering these results in detail, let us examine similar phenomena further 
down the wing span. 

Figure 10 presents pressure distributions at three consecutive mesh intervals centered at the 
y = 0.623 station, and figure 1 1 presents a computed isobar pattern for the upper surface of the 
wing central portion. Every span station used for a mesh plane is identified, thus indicating the 
coarseness of the mesh spacing in the span direction. 

Considering first the sudden pressure jump on the aft portion of the upper surface in figures 9 
and 10. In view of the consistency of this behavior along the entire wing span, it is reasonable to 
regard this as a numerical representation of a rather strong shock. The isobar pattern indicates 
that the discontinuity occurs along a surface inclined about 13° to the free stream. At M - 0.853, 
Cp*(13°) » 0.373; the dashed lines in the figures represent our estimate of the actual shock approx- 
imated by the numerical calculations. The sharp cusp extending into the positive pressure coeffi- 
cient area represents our estimate of a local logarithmic singularity in the streamwise gradient of the 
velocity, which should appear in the solution but which is not well resolved in these calculations 
(ref. 19). 

Consider next the pressure distributions on the upper surface near the wing leading edge in 
figure 10. The curve representing the solution at y/b = 0.623 shows a pronounced dip at x/c « 0.15, 
but the curves at mesh stations on either side show little or no such tendency. Comparisons of suc- 
cessive runs showed that the solution in this region was not stationary; that is, the loop in the isobars 
near the leading edge of station 0.623 shown in figure 1 1 had a spanwise fluctuation with iteration 
count. This can occur in “converged” solutions to difference equations when, for example, some of 
the eigenvalues of the algebraic equations are complex and have absolute values equal to one. 2 We 
hypothesize that this situation has occurred here and attribute it to the coarseness of the mesh. This 
aspect of the interpretation will be raised again in the discussion of theA/oo = 0.825 solution. 


2 The solutions are converged in the sense that further iteration will not significantly reduce the effect on the 
solution of those eigenvalues with modules less than one. 
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The important point to be raised regarding the leading-edge pressure behavior presented in fig- 
ures 9 and 10 is that the abrupt changes at x/c = 0.15 could be attributed to an attempt by the 
numerical procedure to form a shock at this location. At this point, it is impossible to say that a 
converged stationary solution based on a finer mesh would exhibit the same result. However, the 
calculations indicate the strong possibility that wings of the type being studied may have weak 
oblique leading-edge shocks preceding much stronger ones that are nearly normal to the free stream. 


Comparison of Numerical Simulations and Experiments for M ^ = 0.825 

A set of calculated pressure distributions for the upper and lower surfaces of the C-141 wing 
flying at M ^ = 0.825 is shown in figure 12. (The Kutta condition at the trailing edge was satisfied by 
making the upper and lower pressures join at the first point not shown in the figures.) Every fourth 
spanwise station of the moderately fine (200,000 points) mesh calculation is presented. Figure 13(a) 
and (b) show upper surface isobar plots constructed from the coarse and fine mesh computations, 
respectively. For completeness, figure 14 includes the span load distribution obtained from the 
computed circulation for the fine mesh. At two points, the circulation was also calculated from a 
numerical integration of the chord loading. The two methods give the same value within the accuracy 
to which one can estimate the shock location. 

It is appropriate to consider the accuracy of these calculations insofar as they represent a solu- 
tion to equation (32). We can discuss the effect of mesh refinement, the degree of convergence and 
the resolution of the singularity at the shock-wing juncture. 

The general effect of going from an 80,000-point to a 200,000-point mesh is demonstrated by 
comparing (a) and (b) of figure 13. The isobar patterns are significantly smoother for the finer mesh. 
A more detailed visualization is given in figure 1 5, where section pressure distributions are shown for 
both mesh sizes at y/b » 0.4. The better resolution of the peak pressure and shock position for the 
fine mesh calculations is evident. 

The criterion for calling a solution converged in a nonlinear relaxation procedure is not simple. 
It is usually based on the iterative history of certain residuals and is generally considered acceptable 
only when the final result is stationary. In our calculations the deciding parameters were the residuals 
of the span circulation, and we accept the less stringent criterion for convergence given in the preced- 
ing footnote. Under this condition two successive “converged” solutions for the coarse mesh are 
shown in figure 16(a). On the lower surface, both solutions were the same; hence only the upper sur- 
face is shown. Figure 16(b) shows two successive “converged” solutions for the finer mesh. Com- 
pletely stationary convergence is not indicated in either case, although experience with studies of 
two-dimensional flows suggests that the solution found by the finer mesh has settled down to the 
place where the shock location is correctly located to within about one chordwise mesh point, or 
2.5 percent chord for the fine mesh. 

Some final results showing a comparison of flight, wind tunnel, and computed data on the upper 
surface of the C-141 wing are shown in figure 17. The calculations were carried out for an airplane 
flying at 0° angle of attack; that is, for a linearly twisted wing inclined +4° to the free stream at the 
root, and -1° to the free stream at the tip sections. No attempt was made to match flight angle of 
attack, since such a comparison is invalidated by trailing-edge viscous effects. Parameters that 
matched fairly well were leading-edge pressure gradient, peak pressure coefficient, and integrated 
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upper surface loading. The discrepancy in pressure curvature in the region 0.1 < x/c < 0.4 is not 
understood. The diference in shock location given by flight and numerical simulation is within the 
range to be expected on the basis of two-dimensional experience. 

After comparisons were made with the data presented in reference 16, some additional flight 
data on the upper and lower surfaces of the wing and at different span locations were obtained. 
These data are shown in figures 18 and 19 together with the numerical results already discussed. 
Note that the upper surface comparison for M « 0.82 and y/b ~ 0.2 and 0.4 in figure 18(a) and 
(b) are quite good, while the lower surface agreement is quite poor. This is to be expected because 
of the body and nacelle interference on the lower surface at these inboard stations. (It should be 
mentioned that the wing is mounted on top of the fuselage.) The upper surface pressure coefficients 
shown in figure 1 8 indicate that the in-flight shock sweep is greater than that of the computed shock. 
This is understandable because the airplane wing section is varied in the span direction to sweep the 
shock thereby decreasing its strength. In the calculation this span variation of wing section was 
neglected. 


CONCLUSIONS 


Calculations were performed that numerically simulated the three-dimensional, transonic 
( M oo = 0.825) flow field about the C— 141 airplane wing. Computed results compared favorably with 
flight data for cases in which the flow was not separated by embedded shock waves. In fact, the 
computed results predicted the in-flight (Re = 36X10 6 ) pressure distribution more accurately than 
the wind-tunnel tests (Re = 8.5X10 6 ), which showed appreciable trailing-edge separation. Results 
indicate that the inviscid, computational method can be usefully applied to analyze three-dimensional 
transonic flows of engineering interest. The numerical algorithm is currently being extended to treat 
flows with shocks that are more oblique and geometries that include simple wing-body combinations. 


Ames Research Center 

National Aeronautics and Space Administration 
Moffett Field, Calif., 94035, May 21, 1973 
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TABLE 1.- C-141 WING SECTION COORDINATES (STREAMWISE) 

[y/b = 0.389] 


Upper Surface 

Lower Surface 

x/c 

z/c 

x/c 

z/c 

0.0003 

0.0042 

0.0005 

-0.0029 

.0008 

.0059 

.0009 

-.0041 

.0013 

.0073 

.0017 

-.0058 

.0018 

.0084 

.0031 

-.0080 

.0053 

.0131 

.0059 

-.0109 

.0076 

.0154 

.0099 

-.0141 

.0112 

.0182 

.0165 

-.0178 

.0160 

.0213 

.0217 

-.0201 

.0233 

.0252 

.0373 

-.0250 

.0480 

.0343 

.0681 

-.0311 

.0728 

.0408 

.0936 

-.0345 

.0978 

.0458 

.1038 

-.0356 

.1178 

.0492 

.1242 

-.0375 

.1478 

.0534 

.1547 

-.0398 

.1979 

.0589 

.1750 

-.0410 

.2480 

.0630 

.2053 

-.0424 

.2980 

.0659 

.2557 

-.0441 

.3480 

.0678 

.3059 

-.0451 

.3979 

.0688 

.3560 

-.0453 

.4477 

.0688 

.4058 

-.0449 

.4974 

.0678 

.4554 

-.0438 

.5470 

.0659 

.5048 

-.0421 

.5965 

.0630 

.5541 

-.0399 

.6461 

.0591 

.6025 

-.0369 

.6962 

.0543 

.6521 

-.0334 

.7465 

.0485 

.6845 

-.0307 

.7967 

.0414 

.7247 

-.0278 

.8472 

.0329 

.7846 

-.0225 

.8976 

.0235 

.8341 

-.0180 

.9986 

.0017 

.8670 

-.0148 



.8999 

-.0114 



.9164 

-.0096 



.9651 

-.0039 



.9822 

-.0019 
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Figure 1.— Critical pressure coefficient vs. free-stream Mach number given by Eqs. (1) and (14). 
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Figure 5.- Comparison of actual and simulated C-141 wing planform. 
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Figure 6.- Shock-wave pressure jump on airfoils in transonic flow. 
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data for C-141 airplane and model. 


0.85 8.5x10® Tunnel 
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Figure 10.— Pressure coefficient at three successive spanwise mesh stations, = 0.853 coarse mesh. 
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10.— Concluded. 
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(a) Coarse mesh. 

Figure 13.— Isobars for simulation of flow at = 0.825. 
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(b) Fine mesh. 
Figure 13.— Concluded. 
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loading on C-141 wing. 
















Flight data 



(c)y/b «0.63. 
Figure 18.— Concluded. 
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